/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#ifndef __QUANTILE_LOSS_HH__
#define __QUANTILE_LOSS_HH__

#include <algorithm>
#include <boost/optional.hpp>
#include <cmath>
#include <stdint.h>

#include "dyadicp.hh"
#include "frsqrt.hh"
#include "loss.hh"

namespace flassol
{
  class QuantileLoss : public Loss
    {
      private:
        static LossFactory factory;

        float                   tau;
        double                  best_constant;
        double                  best_constant_loss;
        uint64_t                nsamples;
        uint64_t                skipped_update;

      public:
        explicit
        QuantileLoss (float _tau = 0.5f) : tau (_tau),
                                           best_constant (0),
                                           best_constant_loss (0),
                                           nsamples (0),
                                           skipped_update (0)
          { 
          }

        virtual float 
        loss (float p,
              float y)
          {
            best_constant_loss += 
              (y > best_constant) ? tau * (y - best_constant)
                                  : (1.0 - tau) * (best_constant - y);

            ++nsamples;
            float cn = rsqrtf (nsamples);

            best_constant += cn * (y > best_constant ? tau : tau - 1.0);

            return (y > p) ? tau * (y - p) : (1.0 - tau) * (p - y);
          }

        virtual float 
        dldp (float p,
              float y)
          {
            return (y > p) ? -tau : (1.0 - tau);
          }

        virtual std::pair<float, float>
        invariant_update (float p,
                          float wdotx,
                          float adotb,
                          float y,
                          float eta,
                          float anormplusbnorm,
                          float xnorm,
                          bool  adaptive,
                          float adaptivexnorm,
                          float lambda)
          {
            if (y == p)
              {
                return std::make_pair (0.0f, 1.0f);
              }

            float actualeta = (adaptive) ? eta * rsqrtf (xnorm)
                                         : eta / xnorm;
            float xtx = (adaptive) ? adaptivexnorm : xnorm;
            float s = -actualeta * (y > p ? tau : tau - 1.0f);
            float kappa = lambda / (y > p ? tau : tau - 1.0f);

            std::pair<bool, float> ps = 
              flassol::dyadic::p (s,
                                  adotb,
                                  anormplusbnorm,
                                  wdotx,
                                  xtx,
                                  kappa);

            if ((ps.second - y) * (p - y) >= 0.0f)
              {
                // desired step is acceptable

                if (ps.first)
                  {
                    ++skipped_update;

                    return std::make_pair (0.0f, 1.0f);
                  }

                return std::make_pair (s, expf (kappa * s));
              }

            // desired step is too big
            // instead take the maximum acceptable step
            
            std::pair<bool, float> pinv = 
              flassol::dyadic::pinv (y,
                                     adotb,
                                     anormplusbnorm,
                                     wdotx,
                                     xtx,
                                     kappa);

            if (pinv.first)
              {
                ++skipped_update;

                return std::make_pair (0.0f, 1.0f);
              }
            else
              {
                return std::make_pair (pinv.second, expf (kappa * pinv.second));
              }
          }

        virtual std::ostream& 
        report (std::ostream& o,
                LogLevel      level)
          {
            if (nsamples > 0 && level <= INFO)
              {
                o << "best constant (estimated) = " 
                  << best_constant 
                  << std::endl;
                o << "best constant loss (estimated) = " 
                  << (best_constant_loss / nsamples)
                  << std::endl;
              }

            if (skipped_update > 0 && level <= WARN)
              {
                o << "WARNING: skipped " << skipped_update
                  << " updates due to loss of precision "
                  << std::endl;
              }

            return o;
          }
    };
}

#endif // __QUANTILE_LOSS_HH__
